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Abstract 

In this work, we introduce a novel algorithm for the Biot problem based on a Hybrid 
High-Order discretization of the mechanics and a Symmetric Weighted Interior Penalty dis¬ 
cretization of the flow. The method has several assets, including, in particular, the support 
of general polyhedral meshes and arbitrary space approximation order. Our analysis delivers 
stability and error estimates that hold also when the specific storage coefficient vanishes, and 
shows that the constants have only a mild dependence on the heterogeneity of the permeability 
coefficient. Numerical tests demonstrating the performance of the method are provided. 


1 Introduction 


We consider in this work the quasi-static Biot’s consolidation problem describing Darcian flow in 
a deformable saturated porous medium. Our original motivation comes from applications in geo¬ 
sciences, where the support of general polyhedral meshes is crucial, e.g., to handle nonconforming 
interfaces arising from local mesh adaptation or Voronoi elements in the near wellbore region when 
modelling petroleum extraction. Let O c 1 < d ^ 3, denote a bounded connected polyhedral 
domain with boundary dVL and outward normal n. For a given finite time t-p > 0, volumetric load 
f, fluid source g, the Biot problem consists in finding a vector-valued displacement held u and a 
scalar-valued pore pressure held p solution of 

—V-cr{u) + aVp=f in n X (Oitp), (la) 

Codtp + V-(ad(M) — V-(kVp) = g in 17 x (0,7 f)) (lb) 

where cq 0 and a > 0 are real numbers corresponding to the constrained specific storage and 
Biot-Willis coefficients, respectively, k is a real-valued permeability held such that k ^ k a.e. 
in 17 for given real numbers 0 < k ^ k, and the Cauchy stress tensor is given by 

cr{u) := 2/rVsM + Xld^'U, 


with real numbers A 5= 0 and p > 0 corresponding to Lame’s parameters, Vg denoting the sym¬ 
metric part of the gradient operator applied to vector-valued fields, and Id denoting the identity 
matrix of Equations (laI and express, respectively, the mechanical equilibrium and the 

fluid mass balance. We consider, for the sake of simplicity, the following homogeneous boundary 
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conditions: 


M = 0 on (5ri X (Ojtp), (Ic) 

AtVp-n = 0 on (5ri X (0, tp)- (Id) 

Initial conditions are set prescribing «(•, 0) = and, if cq > 0, p{-, 0) = In the incompressible 
case Co = 0, we also need the following compatibility condition on g: 

rg(-,t) = 0 Vte(0,tp), (le) 

Jn 

as well as the following zero-average constraint on p: 

rp(-,t) = 0 Vte(0,tp). (If) 

Jn 

For the derivation of the Biot model we refer to the seminal work of Terzaghi and Biot j^|^. 
A theoretical study of problem Q can be found in |33| . For the precise regularity assumptions 
on the data and on the solution under which our a priori bounds and convergence estimates are 
derived, we refer to Lemma and Theorem |12[ respectively. 

A few simplifications are made to keep the exposition as simple as possible while still retaining all 
the principal difficulties. For the Biot-Willis coefficient we take 


a = 1, 


an assumption often made in practice. For the scalar-valued permeability k, we assume that it 
is piecewise constant on a partition Pq of fl into bounded open polyhedra. The treatment of 
more general permeability coefficients can be done following the ideas of |16| . Also, more general 
boundary conditions than © -( |ld[ ) can be considered up to minor modifications. 

Our focus is here on a novel space discretization for the Biot problem (standard choices are made 
for the time discretization). Several difficulties have to be accounted for in the design of the 
space discretization of problem Q: in the context of nonconforming methods, the linear elasticity 
operator has to be carefully engineered to ensure stability expressed by a discrete counterpart of 
the Korn’s inequality; the Darcy operator has to accomodate rough variations of the permeability 
coefficient; the choice of discrete spaces for the displacement and the pressure must satisfy an inf- 
sup condition to contribute reducing spurious pressure oscillations for small time steps combined 
with small permeabilities when cq = 0. An investigation of the role of the inf-sup condition in the 
context of finite element discretizations can be found, e.g., in Murad and Loula (^[^ . A recent 
work of Rodrigo, Caspar, Hu, and Zikatanov has pointed out that, even for discretization 
methods leading to an inf-sup stable discretization of the Stokes problem in the steady case, 
pressure oscillations can arise owing to a lack of monotonicity. Therein, the authors suggest that 
stabilizing is possible by adding to the mass balance equation an artificial diffusion term with 
coefficient proportional to h^/r (with h and r denoting, respectively, the spatial and temporal 
meshsizes). However, computing the exact amount of stabilization required is in general feasible 
only in 1 space dimension. 


Several space discretization methods for the Biot problem have been considered in the literature. 
Finite element discretizations are discussed, e.g., in the monograph of Lewis and Schrefler |24| ; 
cf. also references therein. A finite volume discretization for the three-dimensional Biot problem 
with discontinuous physical coefficients is considered by Naumovich (^. In 29,30 , Phillips and 


Wheeler propose and analyze an algorithm that models displacements with continuous elements 
and the flow with a mixed method. In j^, the same authors also propose a different method 
where displacements are instead approximated using discontinuous Galerkin methods. In |36| , 
Wheeler, Xue and Yotov study the coupling of multipoint flux discretization for the flow with a 
discontinuous Galerkin discretization of the displacements. While certainly effective on matching 


2 





simplicial meshes, discontinuous Galerkin discretizations of the displacements usually do not allow 
to prove inf-sup stability on general polyhedral meshes. 

In this work, we propose a novel space discretization of problem Q where the linear elasticity 
operator is discretized using the Hybrid High-Order (HHO) method of (c.f. also [I^[T^ 
[Tt]), while the flow relies on the Symmetric Weighted Interior Penalty (SWIP) discontinuous 
Galerkin method of [^, see also [l^ Ghapter 4]. The proposed method has several assets: (i) It 
delivers an inf-sup stable discretization on general meshes including, e.g., polyhedral elements 
and nonmatching interfaces; (ii) it allows to increase the space approximation order to accelerate 
convergence in the presence of (locally) regular solutions; (iii) it is locally conservative on the primal 
mesh, a desirable property for practitioners and key for a posteriori estimates based on equilibrated 
fluxes; (iv) it is robust with respect to the spatial variations of the permeability coefhcient, with 
constants in the error estimates that depend on the square root of the heterogeneity ratio; (v) it 
is (relatively) inexpensive: at the lowest order, after static condensation of element unknowns for 
the displacement, we have 4 (resp. 9) unknowns per face for the displacements + 3 (resp. 4) 
unknowns per element for the pore pressure in 2d (resp. 3d). Finally, the proposed construction 
is valid for arbitrary space dimension, a feature which can be exploited in practice to conceive 
dimension-independent implementations. 

The material is organized as follows. In Section we introduce the discrete setting and formulate 
the method. In Section]^ we derive a priori bounds on the exact solution for regular-in-time vol¬ 
umetric load and mass source. The convergence analysis of the method is carried out in Section]^ 
Implementation details are discussed in Section while numerical tests proposed in Section Fi¬ 
nally, in Appendix]^ we investigate the local conservation properties of the method by identifying 
computable conservative normal tractions and mass fluxes. 


2 Discretization 

In this section we introduce the assumptions on the mesh, define the discrete counterparts of the 
elasticity and Darcy operators and of the hydro-mechanical coupling terms, and formulate the 
discretization method. 


2.1 Mesh and notation 

Denote by H c a countable set of meshsizes having 0 as its unique accumulation point. 
Following [l^ Ghapter 1], we consider h-reflned spatial mesh sequences {Th)heH where, for all 
h e H, Th is a finite collection of nonempty disjoint open polyhedral elements T such that H = 
^ ^ “ m&XTgTy hx with standing for the diameter of the element T. We assume 

that mesh regularity holds in the following sense: For all h e H, Th admits a matching simplicial 
submesh T/i and there exists a real number > 0 independent of h such that, for all h e 7^,(1) for 
all simplex S diameter hs and inradius rg, ghs ^ rs and (ii) for all T sTh, and all S e 

such that S (z T, ghx ^ hs- A mesh sequence with this property is called regular. It is worth 
emphasizing that the simplicial submesh T/i is just an analysis tool, and it is not used in the actual 
construction of the discretization method. These assumptions are essentially analogous to those 
made in the context of other recent methods supporting general meshes; cf., e.g., Section 2.2] 
for the Virtual Element method. For a collection of useful geometric and functional inequalities 
that hold on regular mesh sequences we refer to [T^ Ghapter 1] and . 

Remark 1 (Face degeneration). The above regularity assumptions on the mesh imply that the 
diameter of the mesh faces is uniformly comparable to that of the cell(s) they belong to, i.e., face 
degeneration is not allowed. Face degeneration has been considered, on the other hand, in ^ in the 
context of interior penalty discontinuous Galerkin methods. One could expect that this framework 


3 


could be used herein while adapting accordingly the penalty strategy (131 and (22l. This point lies 
out of the scope of the present work and will be inspected in the future. 


To avoid dealing with jumps of the permeability inside elements, we additionally assume that, for 
all /i e "H, 7ji is compatible with the known partition Pq on which the diffusion coefficient k is 
piecewise constant, so that jumps can only occur at interfaces. 

We define a face P as a hyperplanar closed connected subset of 17 with positive (d—l)-dimensional 
Hausdorff measure and such that(i) either there exist Ti,T 2 e Th such that F c dTi n dT 2 (with 
dTi denoting the boundary of Tf) and F is called an interface or (ii) there exists T e Th such that 
F c dT n dTl and F is called a boundary face. Interfaces are collected in the set boundary 
faces in and we let Fh '■= Fj^yj Fj^. The diameter of a face F e Fh is denoted hy hp- For all 
T e Th, Ft ■= {F e Fh \ F cz dT} denotes the set of faces contained in dT and, for all F e Ft, 
riTF is the unit normal to F pointing out of T. For a regular mesh sequence, the maximum number 
of faces in Ft can be bounded by an integer Ng uniformly in h. For each interface F e F}^, we 
fix once and for all the ordering for the elements Ti,T 2 e Th such that F c dT n dT 2 and we let 
np := riTi.F- For a boundary face, we simply take np = n, the outward unit normal to 17. 

For integers I > 0 and s > 1, we denote by P(;(7/i) the space of fully discontinuous piecewise 
polynomial functions of total degree ^ I on Th and by H^{Th) the space of functions in LfiTl) that 
lie in H^{T) for all T e Th- The notation H^{Pq) will also be used with obvious meaning. Under 
the mesh regularity assumptions detailed above, using [I^ Lemma 1.40] together with the results 
of [19] , one can prove that there exists a real number Capp depending on g and I, but independent 
of h, such that, denoting by tt], the L^-orthogonal projector on ¥}j^{Th), the following holds: For 
all s e {1,...,/+ 1} and all v e F[^{Th), 


< C'app/i'* Vme {0,...,s-1}. (2) 


For an integer / 5= 0, we consider the space 


C\V) -.= c\[o,tF];V), 


spanned by U-valued functions that are I times continuously differentiable in the time inter¬ 
val [Ojtp]. The space C^{V) is a Banach space when equipped with the norm ||:/?||c‘’(y) ■ = 
maxtg|-Q (p] ||:/j(7) IIV, and the space C^{V) is a Banach space when equipped with the norm ||v5||c''(y) ■= 
maxosjmssi l|d™‘/5||co(y)■ For the time discretization, we consider a uniform mesh of the time in¬ 
terval (Ojfp) of step T := tp/N with N e N*, and introduce the discrete times F := nr for 
all 0 ^ n ^ For any if e C^{V), we set fP := if{F) e V, and we introduce the backward 
differencing operator St such that, for all 1 ^ n ^ 


Stf^ := 


If'- — f" 


G V. 


(3) 


In what follows, for X c 17, we respectively denote by (•, ■)j^ and H-Hx the standard inner product 
and norm in T{X), with the convention that the subscript is omitted whenever X = fl. The same 
notation is used in the vector- and tensor-valued cases. For the sake of brevity, throughout the 
paper we will often use the notation a < b for the inequality a ^ Cb with generic constant C > 0 
independent of h, t, cq. A, /x, and k, but possibly depending on g and the polynomial degree k. 
We will name generic constants only in statements or when this helps to follow the proofs. 


2.2 Linear elasticity operator 

The discretization of the linear elasticity operator is based on the Hybrid High-Order method 
of [^. Let a polynomial degree fc > 1 be fixed. The degrees of freedom (DOFs) for the displace- 
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ment are collected in the space 


£/^= X 

iTeTh 


X X v^d-iipy 


(4) 


.FeJ^h 


For a generic collection of DOFs in [7^ we use the notation := (^(vT)TeThy'^F)FeFh) ■ We 
also denote by Vh (not underlined) the function of V^{ThY such that Vh\T = i’t for all T e Th- 
The restrictions of and to an element T are denoted by U^p and = (u^, {vp)p^p^^, 
respectively. For further use, we define the reduction map /J) : 77^(0)'^ ^ [7^ such that, for all 
V e 

i-h'^ = {{FpV)TeThyF%v)p^py), (5) 

where tt^ and tt^ denote the L^-orthogonal projectors on P^(T) and PXi(F'), respectively. For 
all r e 7)i, the reduction map on obtained by a restriction of is denoted by 7^. 

For all r e 77i, we obtain a high-order polynomial reconstruction : [7^ ^ P*+^(T)‘^ of the 
displacement field by solving the following local pure traction problem: For a given local collection 
of DOFs Vt = (vt, (vp)pej^^) e Up, find e such that 

(Vsrp'^^Vp,Vsw)p = (VsVp,Vsw)p + (vp - vp,VsWnpp)p Vw e P^+^(T)‘^. ( 6 ) 


FeJ^t 

In order to uniquely define the solution to ([^, we prescribe the conditions ^p Vp'^^Vp = ^p vp and 


IpVssFp'^^Vp = YiFsFt ^F ^ ® fp — Vp ® npp), where Vgg denotes the skew-symmetric 

part of the gradient operator. We also define the global reconstruction of the displacement : 

Hi 


:d/c+1 


{ThY that, for all G UX, 


K+^Vh)\t = FF^Vt VTgT.. 


(7) 


The following approximation property is proved in 14 Lemma 2]: For all u G i7^(0)'^ni7*+^(PQ)'^, 

||V.(r^+X^u-u)|| (8) 


We next introduce the discrete divergence operator Dp : Up such that, for all q e P^(T) 


{DpVp,q)p = {y■vp,q)p + ^ {vp — vp,qnpp)p (9a) 

FeJ-t 

= —{vp,\^q)p+ ^ {vp,qnpp)p, (9b) 

F^Ft 

where we have used integration by parts to pass to the second line. The divergence operator 
satisfies the following commuting property: For all T e Th and all v G H^{TY, 

D^l^v = ttXV-u). (10) 

The local contribution to the discrete linear elasticity operator is expressed by the bilinear form 
ap on Up X Up such that, for all Wp,Vp G Up, 

aT{Wp,Vp) ■= 2^{{VsFp'^^Wp,Vs1'p'^^'Vp)T + Sp{Wp,Vp)} + \{DpWp,DpVp)T, ( 11 ) 

where the stabilization bilinear form sp is such that 

Spjwp, Vp) := ^ hp^i^ppWp, ^ppVp) p, (12) 

FeFt 
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with face-based residual such that, for all Wj, e U^, 


•= {nFr^^Wrp — wp) — {irpr^^Wrp — wp)- 

The global bilinear form ah on [7* x is assembled element-wise from local contributions: 

ah{Wh,Vh) ■■= Yi aT{wj,,VF). (13) 


Ten 


To account for the zero-displacement boundary condition we consider the subspace 

lllo--={llh={MTen,i^F)FeT,)^UUvF^O G . (14) 


Define on the discrete strain seminorm 


\\v.h\\lh-= Xl W^hWlr, IIH/ille.T := IIVs^’tIIt + X hp^Wvp - vtWf- ( 15 ) 

TeTh FeJ^T 


It can be proved that defines a norm on U_ho- Moreover, using 

the following coercivity and boundedness result for a/,: 
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Corollary 6 ], one has 


V {‘^tJ-)\\Vh\kh ^ \\v.h\\a,h ■= ah{Vh,V.h) ^ ??(2Ai + dA)||:U;^||,_^, 


(16) 


where 77 > 0 is a real number independent of h, r and the physical coefficients. Additionally, we 
know from Theorem 8 ] that, for all w e n H^+‘^{PaY such that V w e H'^+^{Pq) 

and all e U\ g, the following consistency result holds: 

ahiLlw,Vf^) + {V-cT{w),Vh) < {2n\\w\\Hk+2(p^y + A||V-u;||/^fc+i(p„)) \\Vh\\e,h- (17) 

To close this section, we prove the following discrete counterpart of Korn’s inequality. 

Proposition 2 (Discrete Korn’s inequality). There is a real number Ck > 0 depending on g and 
on k but independent of h such that, for all Vf, e g, recalling that Vh e denotes the 

broken polynomial function such that Vh\T = vp for all T e Fh, 

||n/i|| ^ CKdQ\\Vhh,h, (18) 

where da denotes the diameter offl. 

Proof. Using a broken Korn’s inequality on P^{Th)'^ (this is possible since k > 1), one has 

^ -f X + X W^hiFfp^ (19) 


Fen 


Fen 


where Fs,h denotes the broken symmetric gradient on H^{Th)‘^- For an interface F e Fpi 
we have introduced the jump [n/,]p := Thus, using the triangle inequality, we get 

< ||i’F — vtiWf + ll^’F — vt 2 \\f- For a boundary face F e such that F e Ft f Fj^ for 
some T G Th have, on the other hand, ||n/i|p||p = ||np — np||p since vp = 0 (cf. @). Using 
these relations in the right-hand side of (191 and rearranging the sums yields the assertion. □ 


2.3 Darcy operator 

The discretization of the Darcy operator is based on the Symmetric Weighted Interior Penalty 
method of [^, cf. also [I^ Section 4.5]. At each time step, the discrete pore pressure is sought 
in the broken polynomial space 

pk.^^niTh) 
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if Co > 0 , 
if Co = 0 , 


( 20 ) 






where we have introduced the zero-average subspace ■= {qh e ¥’^{Th) \ {qh, 1) = O}. For 

all F G F^, we define the jump and (weighted) average operators such that, for all (p G H^{Th), 
denoting by pT and kt the restrictions of p and k to T e Tk, respectively, 


[p\f '■= PTi — PT2, {y^}F ■= ^TiPTi + ‘^T2^T2, 


( 21 ) 


where ujti = 1 — wt 2 •= 


(kTi +«T2) ' 
2kt-, hT‘ 


Denoting by V/i the broken gradient on H^{Th) and letting. 


for all F G Fk, Ak,_f := ].kt ) ’ define the bilinear form Ch on x Pj^ such that, for all 
qhFh 6 Pu, 


Ch{rh,qh) ■= (kV hTh,'^hqh) — ^ {i{K'Vhrh}F-nF,[qh]F)F + i[rh]F, hqh}F-nF)F) 


FeFi 


Xi ^-l^i['^h]F,[qh]F)F, 


( 22 ) 


FeF: 


Hf 


where c > 0 is a user-defined penalty parameter. The fact that the boundary terms only appear 


on internal faces in (|22|) reflects the Neumann boundary condition (|ld|). From this point on, we 

[13 


will assume that c > with Ctr denoting the constant from the discrete trace inequality 

Eq. (1.37)], which ensures that the bilinear form Ch is coercive (in the numerical tests of Sectionjfi 
we took c = (Ng -I- 0.1)/c^). Since the bilinear form Ch is also symmetric, it defines a seminorm on 
Pjk, denoted hereafter by || j|c,ft (the map || j|c,?i is in fact a norm on 

Remark 3 (Alternative stabilization). To get rid of the dependence of the lower threshold of c, 
on Ctr, one can resort to the BR2 stabilization; c.f. 0 and also Section 5.3.2]. In passing, 
this stabilization could also contribute to handle face degeneration since the penalty parameter no 
longer depends on the inverse of the face diameter (cf. Remark^. This topic will make the object 
of future investigations. 


The following known results will be needed in the analysis. Let 
P* := {r G n H^{Pq) \ n'^r-n = 0 on 5D} , 


P„\ := P^ + Pi;. 


Extending the bilinear form cg to P^y^ x P|^, the following consistency result can be proved 
adapting the arguments of [T^ Chapter 4] to account for the homogeneous Neumann boundary 
condition ([Td|: 

VrGP^,, -{V-{KVr),q) = Ch{r,q) 'iqeP^h- (23) 


Assuming, additionally, that r G P^+^(Pn), as a consequence of 13 Lemma 5.52] together with 


the optimal approximation properties © of TT^ on regular mesh sequences one has. 


sup 

9feeI"d,o('P)\{0} 


Ch{r-TTlr,qh) ^ 

°—■ i’l . ‘ ^ /a 


II II c,h 




(24) 


2.4 Hydro-mechanical coupling 


The hydro-mechanical coupling is realized by means of the bilinear form bh on U\ x P^(7/t) such 
that, for all G U]k and all G F’^iTh), 


bhiFh^qh) ■— ^ bT{vrp,qfi^F): bT{vj.,qh]T) ■——{DFVF,qh\T)T, (25) 

Ten 


where Dif is the discrete divergence operator defined by (9a I. A simple verification shows that, 
for all Vf^ G U’k all qg G P^(77i), 


bh{Vh,qh) ^ \\Vh\\e,h\\qh\\ 


(26) 
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Additionally, using the definition (9a) of and (14) of LLhO': t)® proved that, for all 

- eUi 


jf^ c Q, it holds (xn denotes here the characteristic function of fl), 


bhiVf^,Xn) = 0. 


(27) 


The following inf-sup condition expresses the stability of the hydro-mechanical coupling: 

Lemma 4 (inf-sup condition for bh)- There is a real number /3 depending on fl, g and k but 
independent of h such that, for all e oiTk)) 


llg/ill ^ f3 sup 


bhjVh^qh) 

\\Th\U,h 


(28) 


Proof. Let qh s ^d,oO~h)- Classically 0 . there is e Hq{Q)‘^ such that = q^ and 

\\vq,,\\H^(QY ^ Ik^ill- Let T e Th- Using the iL^-stability of the L^-orthogonal projector (cf., 
e.g., m Corollary 3.7]), it is inferred that 

lIVsTTTUgJlT < ||Vu,J|t. 


Moreover, for all F e Ft, using the boundedness of 7r|l and the continuous trace inequality of 13 
Lemma 1.49] followed by a local Poincare’s inequality for the zero-average function {TT^Vq,, ~'^qh) 
we have 




TTr^qh 


Fh\\F 


S: lIVu 


qh\\T- 


As a result, recalling the definition ([^ of the local reduction map and (15) of the strain norm 
||j|e_T) it follows that < \vq,^\]Ti{T)d. Squaring and summing over T e Fn the latter 

inequality, we get 

ll^ll- (29) 


\\LtVqJe,h 


< IFqhllH^in)'^ ^ 


Using (29), the commuting property (^lOk, and denoting by S the supremum in (28), one has 


khf = i^-Vq^qh) = 2 {DUTVq,,qn)T = -bh{llvq„qn) ^ SWllvq.Un < S||g,.||. □ 

TeTh 


2.5 Formulation of the method 

For all 1 ^ n ^ A^, the discrete solution {uf[,pff) e C/^ o ^ time t" is such that, for all 

{v.h^ qh) e UX.o ^^{'Th), 


ah{ul,Vh) + bh{Vh,pl) = Ihi^-h), (30a) 

{coStPh,qh) - bh{5tul,qh) + Ch{pl,qh) = {g'^,qh), (30b) 

where the linear form l]^ on C/^ is defined as 

ll{vH)--={r,vn)= 2 {r,VT)T. (31) 

T^Th 

In petroleum engineering, the usual way to enforce the initial condition is to compute a displace¬ 
ment from an initial (usually hydrostatic) pressure distribution. For a given scalar-valued initial 
pressure field e L^{Vt), we let ^ := and set with ^ UX o unique solution of 

ahiul,vk = llkk - bhivf„0l) 'iVf, e C/^ q. (32) 


If Co = 0, the value of ^ is only needed to enforce the initial condition on the displacement while, 
if Co > 0, we also set ^ to initialize the discrete pressure. 









Remark 5 (Discrete compatibility condition for cq = 0). Also when cq = 0 it is possible to take 
the test function Qh in (30b I in the full space P^(77i) instead of the zero-average subspace ¥^Q{7h)) 


since the compatibility condition is verified at the discrete level. To check it, it suffices to let 


Qh = Xn in (30b), observe that the right-hand side is equal to zero since 5 " has zero average on Tl 
(cf ([^;, and use the definition \22\ of cu together with ( |27[ ) to prove that the left-hand side also 
vanishes. This remark is crucial to ensure the local conservation properties of the method detailed 
in Section m 


3 Stability analysis 


In this section we study the stability of problem (30) and prove its well-posedness. We recall the 


following discrete Gronwall’s inequality, which is a minor variation of |23[ Lemma 5.1]. 

Lemma 6 (Discrete Gronwall’s inequality). Let an integer N and reals S,G > 0, and K ^ 0 be 
given, and let (a")o-snssAr, {b‘^)oiin^N, and ( 7 ”)osSrassiv denote three sequences of nonnegative real 
numbers such that, for all 0 ^ n N 


+ 6 ^b^ + K ^ 7 ’"a™ + G. 


771 = 0 


771 = 0 


Then, if'y'^6 < 1 for all 0 ^ m ^ N, letting c™ := (1 — 7 ’”( 5 ) it holds, for all 0 ^ n ^ N, 


+ (5 ^ ^ exp (5 ^ c"*!™ x G. 


(33) 


771 = 0 


771 = 0 


Lemma 7 (A priori bounds). Assume f e G ^{L'^(Ll)'^) and g e G^{L^{Ll)), and let (yPh^pl) = 
(iih,p^) with {Uf,,p^) defined as in Section 2.5 For alll ^ n ^ N, denote by {u^,p^) the solution 
to (|30|) . Then, for r small enough, it holds that 


WvLhll.h + Wcfphf + ^ 


N 


\\Ph -Ph ir + 2 ^ ((2^) ^ + Co) b' 


.Oil 2 


Qifj, + dX 

^ 71 = 1 

+ (2^) bQ||/bi(L2(Q)d) + {2fj. + dA)t|||5bo(j;,2(Q)) + Cq ^tFbllco(L2(n))i (34) 
with the convention that Cq = 0 «/co = 0 and, for 0 ^ n N, p^ := [p^, 1). 


Remark 8 (Well-posedness). Owing to linearity, the well-posedness of (30) is an immediate 
consequence of Lemma^ 


Remark 9 (A priori bound for cq = 0). When cq = 0, the choice (20) of the discrete space for 
the pressure ensures that pf, = 0 for all 0 ^ n N. Thus, the third term in the left-hand side 
of (34) yields an estimate on bft ^and the a priori bound reads 


.N\\2 


II 


+ 


1 


N 


2pj -l- d\ 


\\Phf + 1] 


< 


( 2 ^) ^ (c?bl/||ci(L2(n)‘i) + l|P°ll^) + (2^^ + c?A)tp||gbo(i2(Q)). (35) 

The convention Ci)"^ = 0 if cq = 0 is justified since the term T 2 in point (j) of the 

following proof vanishes in this case thanks to the compatibility condition (le). 


Proof of Lemma Throughout the proof, Ci with i e N* will denote a generic positive constant 
independent of h, r, and of the physical parameters cq, A, p,, and k. 
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(1) Estimate of ||pJJ Using the inf-sup condition (281 followed by (271 to infer that 

bhitlhiPh) ~ mechanical equilibrium equation (30aI, and the second inequality in (161, 

for all 1 ^ n ^ -/V we get 

Hh^UIoMO} \\P-h\U,h v^€Ul^\{0} \\V.h\U,h 

= /3 sup ^hiPh)-ah{ul,v^) ^ ^ (2;, + dXy/^ul\\a,H) , 

^TTk \rni ^ 




where we have set, for the sake of brevity, := /3inax(C'K, ?7)- This implies, in particular, 

K ^ 2Ci (4lir IP + (2m + dX)M\\lf,) (36) 


(2) Energy balance. Adding (30a) with (30bI with and summing the 

resulting equation over 1 ^ n ^ A^, it is inferred 


N 


N 


N 


N 


N 


E ra,{ul,S,ul) + ^ ^ r\\pl\\l, = ^ ^ t(5",pD- (37) 


We denote by C and TZ the left- and right-hand side of (371 and proceed to find suitable lower 
and upper bounds, respectively. 

(3) Lower bound for C. Using twice the formula 

2x{x - y) = x'^ + {x - y)'^ - y'^, (38) 

and telescoping out the appropriate summands, the first two terms in the left-hand side of 
can be rewritten as, respectively. 


.n||2 


.0 ||2 


1 1 

=1 n=l “ 

N 11 ^ 1 

X T{coStpl,pl) = -\\cfp^f + b Z1 - ^Wco^plf- 


n=l 

N 


(39) 


n=l 
N\\ 


Using the above relation together with (36l and ||/ || ^ \\f\\c'^{L^{n)‘‘)j it is inferred that 


_ ^ll.,.0||2 I ^\\„^h^N\\2_ i'||„y2„0||2 

- I ^ 


—h \\a,h qII—^ n ITo Ph 


VhW 


N 


1 m2 

_i_ _ ^ _i_ _||„n||2 ^ r \ _ 

^ 8Ci{2p + ^ Kllc,/i- + 


4(2M + dA) MlIcbiUnU) 


)■ (40) 


(4) Upper bound for IZ. For the first term in the right-hand side of (37), discrete integration by 
parts in time yields 


N 


N 


E <{Stui) = (/^,o - if,<) - E r{Str,K-^), 


(41) 


n=l 


n=l 


hence, using the Cauchy-Schwarz inequality, the discrete Korn’s inequality followed by (161 to 
estimate < EMU—^Ha /i for all 1 ^ n ^ (with C 2 '■= C^pf2), and Young’s inequality, one 
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has 


N 




^ g [\mh\\a,h + Wh 


1 ^ 


iiMr'ii^,, 


+ 




AT 


^ ^>"f^r + \\fr + ‘2t^T.r\\sj 


n II2 


(42) 


n=l 


^ ^ iiE + \\ui\\i,h + ^ E 


-h II a,h 


n=0 
cN ||2 I II ^0||2 


^^2^34 ||.|| 

+ - UWcHL^ny)^ 


where we have used the classical bound ||/^^ p + ||/^P + 2tF 2^=1 '^\\^tf'^\\‘^ ^ 

conclude. We proceed to estimate the second term in the right-hand side of (371 by splitting it 

into two contributions as follows (here, g'^ := 1)): 


N 


N 


N 


E r(ff^K) = E r{g\pl-pl) + E <r,Pl) ■■= ^1 + ^2. 


(43) 


n=l 
^11 n||2 


Using the Cauchy-Schwarz inequality, the bound X!”=i '''ll5”P ^ ^F||'?|lco(i 2 (Q)) together with (361 
and Young’s inequality, it is inferred that 


f N r Af 

i^ii< E^ii4ir X E^K- 


1/2 


n=l 


a=l 


1/2 






^f||5||c0(l2(O)) X 1“^ ^{d'hWf^'V + i2p +dX)\\ul\\l^f^)^ 

8 C'it|( 2 /r + dA)||g 4 o(i 2 (a)) + 15(2^ + ^a) 1 ^ Ey 'i' 


(44) 


Owing the compatibility condition (leI, T 2 = 0 if Cq = 0. If Cg > 0, using the Cauchy-Schwarz 
and Young’s inequalities, we have 

I‘l2| < j^F E [ X jtpi E "^llco'Kir [ < y^ll5llco(L=(n)) + ^ E (45) 

[ n=l ) [ n=l ) ^ ^ n = l 


Using (42l, (44), and (45), we infer 


o IKII^ 


N \ . N 2 

,h +E ^iKiiE + ii“°iiE + ^ E ^ii^o'Kir + y^iiffiico(L 

^0 I 2co 


2(n)) 


71=1 


+ 8Ci4(2/r + d\)\\g\\lo^LHn)) + [iQ^2p + dX) 


1 , C'2C'3^ ^2 II f ||2 

-I — j 4ll/llci(L2(n)'*)' (46) 


(5) Conclusion. Using (40), the fact that C = TZ owing to (37), and (46), it is inferred that 

N 


WPih \\a,h + M\<^d^Ph \d + /O,, I n\\ \\Ph “-P(rp + 8 E 


(2// + dA) 




71=1 

N 


p E + § E ^^w^Tpif + G, (47) 

^ 71 = 0 


1 

n=l 
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where Ca '■= inax(l,C'i) while, observing that ||cg'^^p^|| ^ l|co'^^P°|| since is a bounded operator, 
and that it follows from below that + l|p°lP)) 


Cr‘G § (4ll/"lt + l|p”ll“) + 4||c;''-p"lP + 

Z/i Cq 

+ 64C,,J(2,. + + ( 5 ^ 5 ^ + -p 




4II/I 


Ci(L2(0)<J)- 


Using Gronwall’s Lemma jb with a° := and a” := + 4||cg^pJJp for 1 ^ n ^ A^, 


5 := r, 50 := 0 and &" := |p"||2^ ior I ^ n ^ N, K = and 7 " = tk, the 


desired result follows. 


□ 


Proposition 10 (Stability and approximation properties for S°). The initial displacement (321 
satisfies the following stability condition: 


||a°|U..<(2/r)-Vbdo||/°|| + b°ll)- 


(48) 


Additionally, recalling the global reduction map 7^ defined by ([^, and assuming the additional 
regularity po e {Pq), e , and V-m° g i7^+^(Po), it holds 

- Ih“°l|o,/i SI (^‘^h-\\u°\\Hi‘+2(^p^)d + A|| V-M°|bfe + I(pj^) + pPb°lb'= + i(Pr!)) ■ 


Proof. (1) Proof of (48l. Using the first inequality in (I 6 I followed by the definition (32l of uj,, 
we have 


M 


'.Oil 

hWa.h 


< 


sup ...‘ik ab? -- 

Vh^uf o\{o} (2 m) 


= ( 2 a^) 


- 1/2 


sup 


(,„||P|| + iipO,) 


lh\\t,h 


where to conclude we have used the Cauchy-Schwarz and discrete Korn’s ( [18| inequalities for 
the first term in the numerator and the continuity ( [ 2 ^ of by, together with the L^(U)-stability of 
TT^ for the second. (2) Proof of (49l. The proof is analogous to that of point (3) in Lemma 11 


except that we use the approximation properties (|^ of tt^ instead of (54l. For this reason, elliptic 
regularity is not needed. □ 


4 Error analysis 

In this section we carry out the error analysis of the method. 

4.1 Projection 

We consider the error with respect to the sequence of projections {Uh,pp)i^n^N, of the exact 
solution defined as follows: For l^n^N,fPfePl^ solves 

Ch{l^,qh) = Ch{p^,qh) yqh e VdiTh), (50a) 

with the closure condition Once ^ has been computed, e g solves 

ahiul,Vf,) = ll{vy,) - bhivi^,pf,) 'iVf, G t/^ g. (50b) 
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The well-posedness of problems (50a I and ( 50b[ ) follow, respectively, from the coercivity of Ch 
on PJ5 q( 7^) and of ah on C/^g. The projection is chosen so that a convergence rate of 

{k + 1) in space analogous to the one derived in can be proved for the ll-IlQ^^i-norm of the 

displacement at final time t-p. To this purpose, we also need in what follows the following elliptic 
regularity, which holds, e.g., when is convex: There is a real number Ceii > 0 only depending 
on n such that, for all ip e LQ(n), with := {q e \ (g, 1) = O}, the unique function 

C G H^(Vl) n LQ{n) solution of the homogeneous Neumann problem 


is such that 


V-(/cVC) = Ip in n, kSJ( p-n = 0 on dfl, 


(51) 


(52) 


For further insight on the role of the choice (501 and of the elliptic regularity assumption we refer 
to Remark M 

Lemma 11 (Approximation properties for (u^,p^)). Let a time step 1 ^ n ^ iV &e fixed. As¬ 
suming the regularity p" G (Pq), it holds 




(53) 


Moreover, recalling the global reduction map defined by ([^, further assuming the regularity 
m" g , V-m" g H^'^^{Pa), and provided that the elliptic regularity (52) holds, one has 


||^-p"|| (54) 

(2m)V=||m;: - Ilu-\\a,h < (2p||«"||^...(p,). + A||V-«"||p..qp,) + ■ (55) 

with global heterogeneity ratio '■= 

Proof. (1) Proof of (|^. By definition, we have that ||^ - p‘^\\c,h = inf^,^gpfc(p^) \\qh - p‘^\\c,h- To 
prove (53), it suffices to take qt = '^hP'^ right-hand side of the previous expression and use 

the approximation properties (|^ of tt^. 


(2) Proof of (54). Let p G H^{n) solve (51) with ip = p^—p^. From the consistency property (23), 
it follows that 

K - p”ir = -(v-(acVC),^ - p") = Ch{p,p^h - pn = Chip - ^ic,rh - p”). 


Then, using the Cauchy-Schwarz inequality, the estimate (53) together with the approximation 
properties ([^ of tt^, and elliptic regularity, it is inferred that 

wf^h - = Chip - <c,rh -pn^K- ^ickhW^h - pic,h 

< h'=+l7tV^||C||p.(P,)||p"||^..i(P„) < h’^+^p^/^WfPh-p-WmH'^ri^PnP 


and (54) follows. 


(3) Proof of ([5^ . We start by observing that 


SdLh ~ Ah^ 


u 


= sup 


ahjul- Ltu'^,Vh) 

\\ 2 Lh\\a,h 


< sup (56) 

u,,€C/g\{q} (2w ' llH/ilk?! 


where we have used the first inequality in (16). Recalling the definition (311 of the linear form 
the fact that /” = — V-cr(M) + Vp, and using (50a|), it is inferred that 


ahiul- IhU'^,Vh) = ^hiVLh) - a.hilXu'^ ,Vh) - bhiVh^Phn) 

= { - ahiLiu^,Vh) - iW-(Tiu^),Vh)} + {iWp'^,Vh) - bhivh,p^)}. 


(57) 
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Denote by Zi and T 2 the terms in braces. Using it is readily inferred that 

I'Jil ^ (2^\\u^\\f^k+2(^p^y + A||V-M"'||^fc+i(pj^)) (58) 

For the second term, performing an element-wise integration by parts on (Vp, Vh) and recalling 


the definition (251 of and (9aI of Dp with q = p^, it is inferred that 

I'^^l = Xj I ^ -p", (■wf - 1’T)wTi=’)F i 

TeTh I F€Ft j 


(59) 


where the conclusion follows from the Cauchy-Schwarz inequality together with (541. Plug¬ 
ging (58)”(59l into (571 we obtain (55). □ 


4.2 Error equations 


We define the discrete error components as follows: For all 1 ^ n ^ 


p" ._ 

&h ■— “h “/ii 


_ fjn 

Ph ■— Ph Ph- 


(60) 


Owing to the choice of the initial condition detailed in Section 2.5 the inital error (e°,p°) : = 
{Vih ~ ^^Ph ~^h) element in the product space U\ g x P^. On the other hand, for all 

1 ^ n ^ (e)(, p^) solves 


ah{el,v^,) + bu{v^,,pl)=t) yv^eUl 

icoStPh,qh) — bhiSte'^, qn) + Ch{Ph^ qh) = ^hiQh), ^ Ph, 

with consistency error 

^hitth) ■= {g'",qh) - {coStp^,qh) - Ch{^,qh) + bh{StUh,qh)- 


(61a) 

(61b) 

(62) 


4.3 Convergence 


Theorem 12 (Estimate for the discrete errors). Let {u,p) denote the unique solution to (0,/or 
which we assume the regularity 

u^C^{H\Paf)nC\H^+\Pnf), p ^ C\H^+\Pn)) ■ (63) 

If Cq > 0, we further assume p G C^{L^(Lt)). Define, for the sake of brevity, the bounded quantities 
Ml := (2p -I- dX) \\u\\Q2(^pi(^p^y^ + ||cg'^^p||c2(L2(o)d), 


M2 •= 


{2p + dxy/^ 

2p 

+ l|Co'^''p|lco(F'' + i(Ffj))- 


^2/r||tl||ci(F'fe+2(Ps^)d) + A||V-M||ci(pfc+i(pj^)) -I- py^||p||ci(F'= + i(Pri)) 


Then^ assuming the elliptic regularity (52l, it holds, letting \= {p'^,1), 

1 N 


W&h lll/t + \\(^d^Ph ll^ + 2 „ _|_ ~ ^ ^ < {tMi + M~''^M2) ■ (64) 

^ n=l 


Remark 13 (Pressure estimate for cg = 0). In the incompressible case cg = 0, the third term 
in the left-hand side of (64) delivers an estimate on the L^-norm of the pressure since p^ = 0 
(cf ([0;. 
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Proof of Theorem^r^ Throughout the proof, Ci with i e N* will denote a generic positive constant 
independent of h, r, and of the physical parameters cq, A, fi, and k. 

(1) Basic error estimate. Using the inf-sup condition (28l, equation (27l followed by (61a), and 
the second inequality in (16), it is readily seen that 

Wpl-m^P sup = p sup ^ + dA)V^||e-|U,„ 

«heUJo\{2} W^hU.h VheUX„\{0} mhU,h 

(65) 

with = firf^'^. Adding (61a I with to (61b I with = rp^ and summing the 

resulting equation over 1 ^ n ^ A^, it is inferred that 


N 


N 


N 


N 


2 ^ ricM^pl) + ^ t\\pI\1^ = ^ r£l{pl). 


( 66 ) 


Proceeding as in point (3) of the proof of Lemmaj^ and recalling that (e(), p°) = ( 0 , 0), we arrive 
at the following error estimate: 


\'i,h + 


1 


4Ci(2/i-fdA) 


N N 

Ph -Phf + 7)\\cf Phf + T\\plfc,h Z TS'fipD- (67) 


n=l 


n=l 


(2) Bound of the consistency error. Using = codtp" -I- V-(d(M"' — KVp"), the consistency 
property (23), and observing that, using the definition (22) of c/,, integration by parts together 
with the homogeneous displacement boundary condition (|lc|), and ( |2^ , 

Ch{p^-rH,pl) + (v-(d*«"),p)() + hn{5,ul,rh) = 0 , 

we can decompose the right-hand side of (@ as follows: 


N 


N 


N 


Y, rSlipl) = Y ^(co(dtp- - 5trH),Pl) + S rcu{p- - - Tu) 


N 


( 68 ) 


+ T{{yiA,u-),pl-rh) + bh{5trH,pl-rh)}--=Zi+Z2 + Z:,. 


n = l 


For the first term, inserting ±i5tp" into the first factor and using the Cauchy-Schwarz inequality 
followed by the approximation properties of (a consequence of 0) and ( [54| of it is inferred 
that 


N 


1/2 


N 


1/2 


|Ti| < <1 CO r [||d*p" - 5tp-f + ||<5,(p" -^)f] \ X <1 2 rWc^^pl"^ 


n=l 


, n=l 




N 


(69) 


C2 (rVi + /r'=+W2) + - 2 r\\crpU^- 


n=l 


For the second term, the choice (50a) of the pressure projection readily yields 

T 2 = 0. 


(70) 


For the last term, inserting into the first argument of and using the commuting prop¬ 

erty (101 of it is inferred that 

N 


^3 = S ^ S [(V-(dt«- - Stu-),pl - pI)t + - u-),pl - Thh] ■ 

n=l [Ten J 
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Using the Cauchy-Schwarz inequality, the bound — uJ)||t ^ WStiLr'^^ ~ MT)lkT 

valid for all T e Th, and the approximation properties ( |49[ ) and ( |55[ ) of and u]^, respectively, 
we obtain 


N 


1/2 


12:31 s: 2 ^ + \\ StUtu -- uDwl ^] 

l^n=l 

^ C3C1 (rAfi + + 


N 


T,r\\p-k-pU^ 


N 


(71) 


E ll n —‘^ 11 2 

, Ik/i -P/ill ■ 

Using to bound the right-hand side of (|6^, it is inferred 


N 




Ci{2p. + PhW +4 


< 


Ci{2fi + dX) 


N N 

Y,r\\pl-plf + 2Y,r\\c:i^plf + G, (72) 


with G := 4(CiC'3 -I- C 2 ) {rAfi + h^+^Ak)^- The conclusion follows using the discrete Gronwall’s 
inequality @ with S = t, K = \\ef^\\li„ a° = 0 and a” = “ P^V + 2||cf for 

1 ^ n ^ A^, 6" = 4 ||p"|| 2_^, and y’" = 1. □ 


Remark 14 (Role of the choice (50 1 and of elliptic regularity). The choice (50 1 for the projection 
ensures that the term T 2 in step (2) of the proof of Theorem \l^ vanishes. This is a key point to 
obtain an order of convergence of {k + 1) in space. For a different choice, say p^ = this 

term would be of order k, and therefore yield a suboptimal estimate for the terms in the left-hand 
side of (74) below (the estimate (75) would not change and remain optimal). This would also be 
the case if we removed the elliptic regularity assumption (52). 


Remark 15 (BDF2 time discretization). In some of the numerical test cases of Section^ we 
have used a BDF2 time discretization, which corresponds to the backward differencing operator 


);(2), n+2 


n+2 




4(^11 +1 + 

2^ 


(73) 


used in place 0/ ([^. As BDF2 requires two starting values, we perform a first march in time 
using the backward Euler scheme (another possibility would have been to resort to the second- 
order Crank-Nicolson scheme). For the BDF2 time discretization, stability estimates similar to 
those of Lemma can be proved with this initialization, while the error can be shown to scale as 
(compare with The main difference with respect to the present analysis focused 

on the backward Euler scheme is that formula ( |38[ ) is replaced in the proofs by 

2x{3x — iy + z) = x"^ — + (2x — y)^ — {2y — z)^ + (x — 2y -t- z)^. 


The modifications of the proofs are quite classical and are not detailed here for the sake of con¬ 
ciseness (for a pedagogic exposition, one can consult, e.g., l2(Jl\ Chapter 6]). 

Corollary 16 (Convergence). Under the assumptions of Theorem\li^ it holds that 


-f \\cf{p((-p^)\\ + 


1 




N 


1/2 


2 ^\\Ph-P^fc,h 


< 


2/i + dX 

< tNi + h^+^Ak + (’^4) 

tA/i -I- h^^^M2 + h^K^^tp ||p||co(//fc+i(Pj^)). (75) 


, n=l 
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Proof. Using the triangular inequality, recalling the definition (601 of and and (161 of 
ll-lla ^i-norm, it is inferred that 




.k + lrk N 
h i/i“ 


+ (2/i)V=||Vs(r^ij^«^-«^)||, 


Pi: 

\pi 


^\\p^H--pii\\ + \\p^-p^\\ 




< l|Co'>^''ll + l|Co' {Ph 




P^W- 


To conclude, use (641 to estimate the left-most terms in the right-hand sides of the above equations. 
Use (551 and (54), the approximation properties (§ of respectively, for the right-most 

terms. This proves (74). A similar decomposition of the error yields (75). □ 


5 Implementation 

In this section we discuss practical aspects including, in particular, static condensation. The 
implementation is based on the hho platforrrQ which relies on the linear algebra facilities provided 
by the EigenS library [^ . 

The starting point consists in selecting a basis for each of the polynomial spaces appearing in 
the construction. Let s = (si,...,Sd) be a d-dimensional multi-index with the usual notation 
|s|i = X;i=i ® = (^1) Given fc > 0 and T e Th, we denote by a basis 

for the polynomial space (T). In the numerical experiments of Section we have used the set 
of locally scaled monomials: 


: = 


X — Xt 


|s|i ^ k 


(76) 


with Xt denoting the barycenter of T. Similarly, for all F e Fh, we denote by Bp a basis for the 
polynomial space P(5_]^(U) which, in the proposed implementation, is again a set of locally scaled 
monomials similar to (76). 


Remark 17 (Choice of the polynomial bases). The choice of the polynomial bases can have a 
sizeable impact on the conditioning of both the local problems defining the displacement recon¬ 
struction (cf. and the global problem. This is particularly the case when using high 


polynomial orders (typically, 7). The scaled monomial basis (76) is appropriate when dealing 
with isotropic elements. In the presence of anisotropic elements, a better choice is to use for each 
element a local frame aligned with its principal axes of rotation together with normalization factors 
tailored for each direction. A further improvement, originally investigated in in the context of 
dG methods, consists in performing a Gram-Schmidt orthonormalization with respect to a suitably 
selected inner product. In the numerical test cases of Section which focus on isotropic meshes 
and moderate polynomial degrees (k ^Z), the basis (76) proved fully satisfactory. 


Introducing the vector bases Bt ■= {BpY, T e Th, and Bp := {BpY, F e FI,, a basis g for the 


space iZ^,o (HI) is given by 


^h,o ^ ■“ X ■— X ®Fj 

TeTh FeJ^l 


while a basis F^ for the space Pji (cf. @) is obtained setting 


K ■■= X ^T- 

TeTh 


^DL15105 Universite de Montpellier 
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When Co = 0, the zero average constraint in Pj^ can be accounted for using as a Lagrange multiplier 
the characteristic function of Notice also that boundary faces have been excluded from the 
Cartesian product in the definition of to strongly account for boundary conditions. Letting, 
for the sake of brevity, := (^^"), n e N, a simple computation shows that 

dim(W^) = dc&id{Th)N^, dim(W^) = dcaid{Pl)N;^_i, dim{V^) = card(rft)iV^ 

The total DOF count thus yields 

dcaTd{Th)N^ + dcaTd{Pl)N^_^ + card(r^)iV^ (77) 


In what follows, for a given time step 0 ^ n ^ N, we denote by and the vectors collecting 
element-based and face-based displacement DOFs, respectively, and by P” the vector collecting 
pressure DOFs. 


Denote now by A and B, respectively, the matrices that represent the bilinear forms Uh (cf. 
and bh (cf. (25l) in the selected basis. Distinguishing element-based and face-based displacement 
DOFs, the matrices A and B display the following block structure: 


bXTd.^J.P. 


B = 


Br 

B.F 


For every mesh element T e Th, the element-based displacement DOFs are only coupled with 
those face-based displacement DOFs that lie on the boundary of T and with the (element-based) 
pressure DOFs in T. This translates into the fact that the submatrix A 7 - 7 - is block-diagonal, i.e.. 


Arr = diag(ATr)Terhi 

with each elementary block of size dim(By)^. Additionally, it can be proved that the blocks 
kxT, T e Th, are invertible, so that the inverse of A 7 - 7 - can be efficiently computed setting 

A7-7- = diag(Aj..j,) 7 ’eTh- (f8) 


The above remark can be exploited in practice to efficiently eliminate the element-based displace¬ 
ment DOFs from the global system. This process, usually referred to as “static condensation”, is 
detailed in what follows. 


For a given time step 1 ^ n ^ A^, the linear system corresponding to the discrete problem (301 is 
of the form 


Arr:Arr; Br 

Arr Arr Br 




pn 

Or 

(79) 

-B^ -Bj. ^C-fcoM 


pn 


G” 



where C denotes the matrix that represents the bilinear form Ch in the selected basis, M is the 
(block diagonal) pressure mass matrix, FVj- is the vector corresponding to the discretization of the 
volumetric load while is the zero vector of length dim(l/^). Denoting by G" the vector 
corresponding to the discretization of the fluid source g”, when the backward Euler method is 
used to march in time, we let 9 = 1 and set 

G" := rG" - BU”"^ 


For the BDF2 method (and n > 2) , we let 0 = 3/2 and set 

G" := ^rG” - ^BU""^ - ^BU""^ 
3 3 3 
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Recalling (78), instead of assemblying the full system, we can effectively compute the Schur com¬ 
plement of A 7 - 7 - and code, instead, the following reduced version, where the element-based dis¬ 
placement DOFs collected in the subvector \J^ no longer appear: 




B_f A^jtA^^Bt- 
jC -f- cqM -|- B^A-^-^^Bt" 



ru^i 



pn 



— A"'" A~^ F" 

G" + b:^ A:^^F^ 


(80) 


All matrix products appearing in ( |M| ) are directly assembled from their local counterparts (i.e., the 
factors need not be constructed separately). Specifically, introducing, for all T e TTt, the following 
local matrices A(T) and B(T) representing the local bilinear forms ut (cf. ( |TT| )) and bx (cf. ([^), 
respectively: 

Att : ^tFt 


A{T) = 






TeTh 






TeTh 


BxAj,}pATj^T , 


and, similarly, for the right-hand side vector 


aT a —1 
^TT^TT^T 


TeTh 


AT a— 1pn 

'^TTt^TT'^Tj 


B^ 


one has for the left-hand side matrix, denoting by 
a global DOF map. 


TeTh 


B(T) = 

the usual assembly procedure based on 


Arj^'Arr^T 


TeTh 


B^A:^^ 


Br 


TeTh 


At ^tt ® > 

ByAjiy Bt, 




TeTh 




The advantage of implementing (80) over (79) is that the number of DOFs appearing in the linear 
system reduces to (compare with (77)) 


dcavA{F),)N^_^ + cwd{Th)N^. 


(81) 


Additionally, since the reduced left-hand side matrix in ( |80[ ) does not depend on the time step n, 
it can be assembled (and, possibly, factored) once and for all in a preliminary stage, thus leading 
to a further reduction in the computational cost. Finally, for all T e Th, the local vector of 
element-based displacement DOFs can be recovered from the local right-hand side vector FJ and 
the local vector of face-based displacement DOFs and (element-based) pressure DOFs (U^^,P^) 
by the following element-by-element post-processing: 


UJ = Ay^ (FJ - At^^U^^ - BtP?^) . 


6 Numerical tests 

In this section we present a comprehensive set of numerical tests to assess the properties of our 
method. 


6.1 Convergence 


We first consider a manufactured regular exact solution to confirm the convergence rates predicted 
in (64). Specifically, we solve the two-dimensional incompressible Biot problem (cq = 0) in the 
unit square domain Q, = (0,1)^ with tp = 1 and physical parameters /r = 1, A = 1, and n = 1. 
The exact displacement u and exact pressure p are given by, respectively 


u{x,t) = ( — sin(7rt) cos(7rxi) cos(7ra:2), sin(7rt) sin(7ra:i) sin(7rT2)), 
p{x,t) = — cos(7rt) sin(7ra:i) cos(7ra;2)- 
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Figure 1: Triangular, hexagonal-dominant, Voronoi, and nonmatching quadrangular meshes for the numerical tests. 
The triangular and nonmatching quadrangular meshes were originally proposed for the FVCA5 benchmark [22] , 
whereas the hexagonal-dominant mesh is the same used in |18| Section 4.2.3]. 


The volumetric load is given by 

f{x,t) = 67r^(sin(7ri) + 7rcos(7rt)) x ( — cos(7ra;i) cos(7ra;2), sin(7ra;i) sin(7ra;2)), 

while g{x,t) = 0. Dirichlet boundary conditions for the displacement and Neumann boundary 
conditions for the pressure are inferred from exact solutions to dil. 


We consider the triangular, (predominantly) hexagonal, Voronoi, and nonmatching quadrangular 
mesh families depicted in Figure The Voronoi mesh family was obtained using the PolyMesher 
algorithm of |34| . The nonmatching mesh is simply meant to show that the method supports 
nonconforming interfaces: refining in the corner has no particular meaning for the selected solution. 
The time discretization is based on the second order Backward Differentiation Formula (BDF2); 
cf. Remark 15 The time step r on the coarsest mesh is taken to be 0.1/2^^ for every choice of 


the spatial degree k, and it decreases with the mesh size h according to the theoretical convergence 
rates, thus, if /12 = hi/2, then T 2 = ri/2^T^. Figure [^displays convergence results for the various 
mesh families and polynomial degrees up to 3. The error measures are \\Ph pressure 

and — I^u^\a,h for the displacement. Using the triangle inequality together with (64) and 
the approximation properties 0 of TT^ and of °Lh)j it i® ^ simple matter to prove that 

these quantities have the same convergence behaviour as the terms in the left-hand side of (64). 
In all the cases, the numerical results show asymptotic convergence rates that are in agreement 
with theoretical predictions. This test was also used to numerically check that the mechanical 
equilibrium and mass conservation relations of Lemma |18| hold up to machine precision. 


The convergence in time was also separately checked considering the method with spatial degree 
A: = 3 on the hexagonal mesh with mesh size h = 0.0172 and time step decreasing from r = 0.1 to 
T = 0.0125. With this choice, the time-component of the error is dominant, and Figure [^confirms 
the second order convergence of the BDF2 scheme. 


6.2 Barry and Mercer’s test case 

A test case more representative of actual physical configurations is that of Barry and Mercer [^, 
for which an exact solution is available (we refer to the cited paper and also to Section 4.2.1] 
for the its expression). We let U = (0,1)^ and consider the following time-independent boundary 
conditions on dU 

UT = 0, ri^Vun = 0, P = 0, 

where r denotes the tangent vector on dfl. The evolution of the displacement and pressure fields 
is driven by a periodic pointwise source (mimicking a well) located at Xq = (0.25,0.25): 

g = S{x — Xq) sin(t), 

with normalized time t := fit for /3 := (A + 2/i)k. As in 
physical parameters: 

co = 0, F;=l•10^ v = 0.1, «: = l•10"^ 


30 


32 


, we use the following values for the 
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Figure 3: Time convergence rate with BDF2, hexagonal mesh 



(a) i = ■^/2 
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(b) i = 31^2 
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jle+5 
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iO 
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:-le+5 
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Figure 4: Pressure field on the deformed domain at different times for the finest Cartesian mesh containing 4,192 
elements 


where E and v denote Young’s modulus and Poisson ratio, respectively, and 

El/ E 

^ (l + i/)(l-2r.)’ 2(1 + 1/) ■ 


In the injection phase i e ( 0, tt), we observe an inflation of the domain, which reaches its maximum 
at i = cf. Figure 


4a 


In the extraction phase i e (tt, 27r), on the other hand, we have a 


contraction of the domain which reaches its maximum at i = 3-n-/2; cf. Figure 


4b 


The following results have been obtained with the lowest-order version of the method corresponding 
to fc = 1 (taking advantage of higher orders would require local mesh refinement, which is out of 
the scope of the present work). In Figure]^ we plot the pressure profile at normalized times i = ’^2 
and t = 3/7'2 along the diagonal (0,0)-(l,l) of the domain. We consider two Cartesian meshes 
containing 1,024 and 4,096 elements, respectively, as well as two (predominantly) hexagonal meshes 
containing 1,072 and 4,192 elements, respectively. In all the cases, a time step r = (^’t'/s) • 10“^ is 
used. We note that the behaviour of the pressure is well-captured even on the coarsest meshes. 
For the finest hexagonal mesh, the relative error on the pressure in the L^-norm at times t = 2 

and t = 3^2 is 2.85%. 


To check the robustness of the method with respect to pressure oscillations for small permeabilities 
combined with small time steps, we also show in Figure]^ the pressure profile after one and two 
step with K = 1 ■ 10“® and r = 1 • 10“^ on the Cartesian and hexagonal meshes with 4,096 and 
4,192 elements, respectively. We remark that the first time step is performed using the backward 
Euler scheme, while the second with the second order BDF2 scheme. This situation corresponds 
to the one considered in Figure 5.10] to highlight the onset of spurious oscillations in the 
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Figure 5: Pressure profiles along the diagonal (0, 0)—(1,1) of the domain for different normalized times t and meshes 
{k ^ 1). The time step is here r ^ ■ 10“^. 


pressure. In our case, small oscillations can be observed for the Cartesian mesh (cf. Figure [fe 
and Figure 6cI, whereas no sign of oscillations in present for the hexagonal mesh (cf. Figure phr 
and Figure 6dl. One possible conjecture is that increasing the number of element faces contributes 
to the monotonicity of the scheme. 


Acknowledgements The work of M. Botti was partially supported by Labex NUMEV (ANR- 
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A Flux formulation 


In this section, we reformulate the discrete problem (301 to unveil the local conservation properties 
of the method. Before doing so, we need to introduce a few operators and some notation to treat 
the boundary terms. 


We start from the mechanical equilibrium. Let an element T e 77i be fixed and denote by Ust '■ = 
the broken polynomial space of degree ^ A: on the boundary dT of T. We define the 
boundary operator : Ugx Ugx such that, for all ip e UgT, 


Lxip\F ■= '^F (v’lF - {ip\p)FeFT)) ^ Tt- (82) 


We also need the adjoint 


of Lp such that 


yipeUgr, {L^Tip,i^)gT = {‘P,L^T*i’)eT (83) 


For a collection of DOFs Vp e U_ti we denote in what follows by vgp e Ugr the function in Ust 
such that vgT\F = vf for all F e Ft- Finally, it is convenient to define the discrete stress operator 


StUt sF^^Vjf + XldD^Vp. 


(84) 


To reformulate the mass conservation equation, we need to introduce the classical lifting operator 
h ■ ^ such that, for all qgB Pj^, it holds 

(.Rihqh,^H)= T. i[Qh]F,{<h}F-nF)F V^eP^Hr.)'". (85) 
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(a) Cartesian mesh (card(7^) = 4,028), first step 
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(b) Hexagonal mesh (card(7)i) = 4,192), first step 
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(c) Cartesian mesh (card(7)^) = 4,028), second step 
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(d) Hexagonal mesh (card(7h,) = 4,192), second step 


Figure 6: Pressure profiles along the diagonal (0, 0)-(l, 1) of the domain for = 1 • 10~® and time step r = 1 • 10““^. 
Small oscillations are present on the Cartesian mesh (left), whereas no sign of oscillations is present on the hexagonal 
mesh (right). 
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Lemma 18 (Flux formulation of problem (301). Problem (301 can be reformulated as follows: 
Find {ufhTP'h) £ Uji 0 ^ smc/i that it holds, for all Qh) G UX o ^ ^dO~h) all T e 7h, 

- Phid, VsVt)t + 2 (^TF(MT'iKir)j I’-F - vt)f = (/", vt)t, 
FeJ-t 

( 86 a) 

{coStPh,qh)T - {StUr - K{\7hPl - Ri^hPDFhqh)T - {4>TF{dtuf-,pl),qh\T)F = {g'",qh), 


FeFt 


( 86 b) 


where, for all T elf and all F e J-t, the numerical traction ^^f ■ H-T ^ ^diF 
mass flux (j)^p : V’^_i{F)‘^ x P^(7)i) ^ Pd-i('^) 

^TFiP-TJ q) •” {^tH-T ~ qId)f^TF + (2/i)-Lj,’ (f)gpLp{vgT ~ 


4‘TfFfi qh) '■ — 


_f + {KV;ig;i}F)-'n.TF — ^^^[9 ;i]f(»T’TF-«f) if F e Tf, (87) 


otherwise. 


with t) 5 T G IPd(-^7’) SMc/i that \]dT\F = ^f for all F e J-'t, ond it holds, for all F e Tf such that 
F e Rti gi Ft 2 , 


^Tif(Mt'i)P/i|Ti) + ^T2F(l*T’2’7'/t|T2) “ 0 
^TiFidtUp,pf) + cj)p^p{6tUp,pf) = 0. 


( 88 a) 

( 88 b) 


Proof. (1) Proof of ( 86 al. Proceeding as in [10| Section 3.1], the stabilization bilinear form st 


defined by (121 can be rewritten as 

St{Wtj2Lt) = Xj i^T*i^dT^Ti'^dT — Wt)),Vp — Vt)f- 


FbFt 


Therefore, using the definitions ([^ of r^'^^Vp with w = and (9al of DpVp with q = pf\T, 

and recalling the definition ([M| oi Sp, one has 


aT{yfj-,Vp) = (S'^M^,Vsr’T)r+X (■S'tMt’^tf + (2^)Ly*(f)5rLT(“5T~“T'))i '*^f —'1’t)f- (89) 

F^Ft 


(90) 


On the other hand, using again the definition ( |9a| ) of DpVp with q = pf\Ti one has 

bT{Vp,p^\T) = -{PhIdFsVT)T - X (PMT^TFj’WF -'1’t)f- 

FeFt 


Equation ( 86 a I follows summing (891 and (90). 

(2) Proof of ( 86 b). Using the definition (9b I of Dp with Vp = StUp and q = qh\T! it is inferred 
that 

briStUT^qh) = —idtUp,'\/hqh)T + X {dtvff-nTF,qh\T)F- (91) 

FeFt 

On the other hand, adapting the results [l^ Section 4.5.5] to the homogeneous Neumann boundary 
condition (jld)), it is inferred 


Ch{pl,qh)= X iFFhPl-RidPiyyhqhh 


T^Th 


X ({s^V^tP/dF-^TF - ^-^^^^[K]f(wtF-Wf), g/i|T)F i- (92) 


FeFtf\F\ 
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Equation (86b I follows summing (91) and (92). 


(3) Proof of (88). To prove (88aI, let an internal face F e be fixed, and make in (88a) such 
that vt = 0 for all T e Th, vpi = 0 for all F' e Fh\{F}, let vj^ span P^_]^(F) and rearrange the 
sums. The mass flux conservation (88b) follows immediately from the expression of (j)!fp observing 


that, for all qh) e U_l x Pj^ and all F B the quantity 


[ — vp + {kV hqh}F)-nF — [Qh\F 


is single-valued on F. 


□ 


Let now an element T e Th he fixed. Choosing as test functions in ( |86a| ) e Uff g such that 
vp = 0 for all F e Th, vp' = 0 for all T' e Th\{T}, and vp spans P^(T)‘^, we infer the following 
local mechanical equilibrium relation: For all Vp e P^(T)‘^, 

iS'fyJf-plIh,V,vp)p- ^ i^^ppiuJf,pl\p),vp)p = {r,vp)p. 


Similarly, selecting qh in (86b) such that qh\p' = 0 for all T' e Th\{T} and qp := qh\p spans P^(T), 
we infer the following local mass conservation relation: For all qp e P^(r), 

{co5tpl,qT)T - {5tU^ - k{V hPl - R’l,hPh)F<lT)T - XI {4>PF{^tUfr,pD,qT)F = {g'^,qT)- 

FeFt 


To actually compute the numerical fluxes defined by (87), besides the operator Sp defined by 
(which is readily available once and Dp have been computed; cf. ([^ and ( [9a| , respectively), 
one also needs to compute the operators Lp and i^’*. The latter operation can be performed at 
marginal cost, since it only requires to invert the face mass matrices of P^_j^(F) for all F e Tp. 
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